#include "dEdxGausN.h"

Double_t Wproton   (){ return  0.075;}
Double_t Wpion     (){ return  0.075;}
Double_t Wkaon     (){ return  0.075;}
Double_t Welectron (){ return  0.075;}
Double_t Wmuon     (){ return  0.075;}

//toy x=10 GeV/c from mn2d.C 10e test
Int_t Nproton(){   return 30000;}
Int_t Npion(){     return 49000;}
Int_t Nkaon(){     return 20000;}
Int_t Nelectron(){ return  1000;}
Int_t Nmuon(){     return   500;}

TRandom3 fran(1);

void FillHist(TH2D * hh, const Int_t ibin, const Double_t nsample, const Double_t mass, const Double_t ww)
{

  Double_t bbpar[10];
  memcpy(bbpar, dEdxGausN::IniBBpar(), sizeof(Double_t)* dEdxGausN::NparBB());

  const Double_t xx = hh->GetBinCenter(ibin);
  const Double_t my = dEdxGausN::GetMean(xx, mass, bbpar);

  for(Int_t ip=0; ip<nsample; ip++){
    hh->Fill(xx,fran.Gaus(my, my*ww));
  }
}

void genhist()
{
  TList *ll=new TList;

  TH2D * horiginal = dEdxGausN::GetHistTPCsigVsLog10P("horiginal");
  ll->Add(horiginal);
 
  for(Int_t ibin = 1; ibin<=horiginal->GetNbinsX(); ibin++){
    FillHist(horiginal, ibin, Nproton(),   dEdxGausN::Mproton(),   Wproton());
    
    //const Double_t xx = horiginal->GetBinCenter(ibin);
    //-1 ~ 2 -> 0~pi
    //const Double_t arg = (xx+1)*6;
    //TF1 *f1=new TF1("f1","49*(1+0.5*TMath::Sin((x+1)*6))/(49*(1+0.5*TMath::Sin((x+1)*6))+30+20+1)",-1,2)
    //FillHist(horiginal, ibin, Npion()*(1+0.5*TMath::Sin(arg)),     dEdxGausN::Mpion(),     Wpion());

    FillHist(horiginal, ibin, Npion(),     dEdxGausN::Mpion(),     Wpion());
    
    if(dEdxGausN::Ntype()>2) 
    FillHist(horiginal, ibin, Nkaon(),     dEdxGausN::Mkaon(),     Wkaon());
    
    if(dEdxGausN::Ntype()>3) 
    FillHist(horiginal, ibin, Nelectron(), dEdxGausN::Melectron(), Welectron());

    if(dEdxGausN::Ntype()>4) 
    FillHist(horiginal, ibin, Nmuon(),     dEdxGausN::Mmuon(), Wmuon());
  }

  TH1D *horiginalyield = horiginal->ProjectionX(Form("%syield", horiginal->GetName()), 1, horiginal->GetNbinsY(),"oe");
  ll->Add(horiginalyield);

  TH2D *hreduced = (TH2D*) horiginal->Clone("hreduced");
  ll->Add(hreduced);
 
  TH1D *hreducedyield = hreduced->ProjectionX(Form("%syield", hreduced->GetName()), 1, hreduced->GetNbinsY(),"oe");
  ll->Add(hreducedyield);

  TFile *fout = new TFile(Form("outgenhist_TOY%d.root", dEdxGausN::Ntype()),"recreate");
  ll->Write();
  fout->Save();
  fout->Close();
}

int main(int argc, char * argv[])
{
  for(int ii=0; ii<argc; ii++){
    printf("%d: %s\n", ii, argv[ii]);
  }
  if(argc!=1){
    printf("argc!=1\n");return 1;
  }

  dEdxGausN::Ini();

  genhist();
  return 0;
}
